#ifndef UTILS_HPP
#define UTILS_HPP
#include <armadillo>
#include <vector>

static const double PI=3.1415926535;

using namespace arma;
using namespace std;


template<typename T>
mat vect2mat( T *vect,int n, int p)
{
	mat A(n,p);
	int k=0;
	for(int i=0;i<p;i++)
	{
		for(int j=0;j<n;j++)
		{
			A(j,i)=(double)vect[k];
			k++;
		}
	}
return A;	
}
template<typename T,typename T2>
T2 vect2mat( T *vect,int n, int p)
{
	T2 A(n,p);
	int k=0;
	for(int i=0;i<p;i++)
	{
		for(int j=0;j<n;j++)
		{
			A(j,i)=(T)vect[k];
			k++;
		}
	}
return A;	
}
	
template<typename T>
void  mat2vect( mat A,int n, int p, T* vect)
{
	int k=0;
	for(int i=0;i<p;i++)
	{
		for(int j=0;j<n;j++)
		{
			vect[k]=(T)A(j,i);
			k++;
		}
	}

}
	
template<typename T>
void  spmat2vect( mat A,int n, int p, T* vect)
{
	int k=0;
	for(int i=0;i<p;i++)
	{
		for(int j=0;j<n;j++)
		{
			vect[k]=(T)A(j,i);
			k++;
		}
	}

}

mat growingvect(int n)
{
	mat X(n,1);
	for(int i=0;i<n;i++)
	{
		X(i,0)=i;
	}
	return X;
}

void Arangemat(mat &X, mat v)
{
	int m=X.n_rows;
	int n=X.n_cols;
	mat temp(m,n);
	for(int i=0;i<m;i++)
	{
		temp(i,span::all)=X(v(i,0),span::all);
	}
	X=temp;		
}
void Sort(double *w,int N, mat &X)//Bubble sort
{
    double y=0;
    double foo=0;
    bool swapped = true;
    while(swapped)
    {
        swapped = false;
        for(int i=1;i<N;i++)
        {
            if(w[i-1]>w[i])
            {
    		y=w[i];
                w[i]=w[i-1];
                w[i-1]=y;
                foo=X(i,0);
                X(i,0)=X(i-1,0);
                X(i-1,0)=foo;
                swapped=true;
            }
        }

    }

}
template<typename T> 
void Sort( T &w)//Bubble sort
{
    int y=0;
    bool swapped = true;
    int N=w.n_elem;
    while(swapped)
    {
        swapped = false;
        for(int i=1;i<N;i++)
        {
            if(w(i-1)>w(i))
            {
    		y=w(i);
                w(i)=w(i-1);
                w(i-1)=y;
                swapped=true;
            }
        }

    }

}



template<typename T>
T min2(T a,T b)
{
	if(a<b)
	{
		return a;
	}else{
		return b;
	}
}

//Binary Search
//


template<typename T>
T abs2(T a)
{
	if(a>0)
	{
		return a;
	}else{
		return -a;
	}
}
inline double Phi(double x)
{
	boost::math::normal_distribution<> d(0,1); 
	double t=cdf(d,x);
	if(t==1)
	{
		t=0.99999;
	}
	if(t==0)
	{
		t=0.00001;
	}
	return t;
	//return boost::math::cdf<double, boost::math::policies::policy<> >(boost::math::normal_distribution(0,1),x);
}

mat add_mat(mat m1 ,mat m2,int p)
{
	mat m3(p,1);
	m3.zeros();
	int k=0;
	for(int i=0;i<p;i++)
	{
		if(m2(i,0)!=0)
		{
			m3(i,0)=m1(k,0);
			k++;				
		}	
	}
	return m3;
}

	mat vectortomat(vector<mat> Pi,int m)
	{
		int sum_n=0;
		for(int i=0;i<m;i++)
		{
			sum_n+=Pi[i].n_rows;

		}
		int p=Pi[0].n_cols;
		mat theta(sum_n,p);

		int sumb=0;
		for(int i=0;i<m;i++)
		{
			if(Pi[i].n_rows!=0)
			{
			theta(span(sumb,sumb+Pi[i].n_rows-1),span::all)=Pi[i];
			sumb+=Pi[i].n_rows;
			}
		}
		return theta;

	}


 
void   vec2vec(vector<double*> Pi,int m,int n,double *theta)
{

		int k=0;
		for(int i=0;i<m;i++)
		{
			for(int j=0;j<n;j++)
			{

				//cout << "\n" << Pi[i][j];
				theta[k]=Pi[i][j];
				k++;			
			}	
		}
}

template<typename T>
void mat2vector(mat X, T* Y)
{

	int n=X.n_rows;
	for(int i=0;i<n;i++)
	{
		Y[i]=X(i,0);	
	}
}

mat add_one(mat X,int a)
{
	int n=X.n_cols;
	int m=X.n_rows;
	mat X2(X.n_rows,n+1);
	for(int i=0;i<m;i++)
	{
		X2(i,n)=1;
	}
	X2(span::all,span(0,n-1))=X;
	return X2;
}

mat  del_one(mat X, int a)
{
	int n=X.n_rows;
	int p=X.n_cols;
	mat Y(n,p-1);
	int k=0;
	for(int j=0;j<p;j++)
	{
		if(j!=a)
		{
			Y(span::all,k)=X(span::all,j);
			k++;
		}
		
	}
	return Y;
}

template<typename V>
int is_present(V X,double t)
{
	int n=X.n_elem;
	for(int i=0;i<n;i++)
	{
		if(X(i)==t)
		{
			return 1;
		}
	}
	return 0;
}

template<typename T>
T merge(T X1,T X2)
{
	int n1=X1.n_elem;	
	int n2=X2.n_elem;
	int n=n1+n2;	
	if(n!=0)
	{
		T X(n);
		int k=0;
		for(int i=0;i<n1;i++)
		{
			X(k)=X1(i);
			k++;
		}
		for(int i=0;i<n2;i++)
		{
			X(k)=X2(i);
			k++;
		}
		Sort<T>(X);
		return X;
	}else{
		return X1;
	}
}

template<typename T>
T rbind(T X1, T X2)
{
	int p1=X1.n_cols;
	int p2=X2.n_cols;
	int n1=X1.n_rows;
	int n2=X2.n_rows;
	if(p1!=p2){ cout << "Warning p1!=p2";}
	T X(n1+n2,p1);
	int k=0;
	for(int i=0;i<n1;i++)
	{
		X(k,span::all)=X1(i,span::all);
		k++;
	}
	for(int i=0;i<n2;i++)
	{
		X(k,span::all)=X2(i,span::all);		
		k++;
	}
	return X;


}
template<typename T>
void vect2vect(vector<T> X, T* res)
{
	int n=X.size();
	for(int i=0;i<n;i++)
	{
		res[i]=X[i];
	}
}


mat del_lines(mat X, int index)
{
	int n=X.n_rows;
	int p=X.n_cols;
	mat f=(sum(X,0));
	int max=f(0,index);
	mat Y(max,p);
	int k=0;
	for(int i=0;i<n;i++)
	{
		if(X(i,index)==1)
		{
			Y(k,span::all)=X(i,span::all);
			k++;
		}
	}
	cout << Y;
	return Y;
}
int max_ones(mat X,double& max)
{
	int p=X.n_cols;
	mat foo=mean(X,0);
	int i=0;
	int k=0;

	while(foo(0,i)==1)
	{
		i++;
	}
	k=i;
	double swap=foo(0,i);	

	for( int j=i+1;j<p;j++)
	{
		if((foo(0,j)>swap) & (foo(0,j)!=1))
		{
			k=j;
			swap=foo(0,j);
		}
	}
	max=swap;
return k;

}

vector<vector<int> > GreedySearch(mat X,double thres)
{
	int T=1;
	vector<int> resp;
	vector<int> resb;
	double max=0;
	mat Y=X;
	int p=X.n_cols;
	int index=0;
	mat foo=mean(X,0);

	for(int i=0;i<p;i++)
	{
		if(foo(0,i)==1)
		{
			resp.push_back(i);
			resb.push_back(1);

		}	
	}
	while(T)
	{
		mat X1=Y;
		index=max_ones(X1,max);//returns max mean columns other than one// max is the value of this max
		cout << index;
		cout << X1;
		Y=del_lines(X1,index);
		if(max<thres)
		{
			T=0;
		}else{
			resp.push_back(index);
			resb.push_back(max);
		}
	}
	vector<vector<int> > res;
	res.push_back(resp);
	res.push_back(resb);
return res;
}

template <typename T>
double Norm(T X,T Y)
{
	double sum;
	int p=X.n_elem;
	for(int j=0;j<p;j++)
	{
		sum+=(double)pow(X(j)-Y(j),2);
	}
	return sum;
}
double logfact(int i)
{
	double res=0;
	while(i>1)
	{
		res+=log(i);
		i--;
	}
	return res;
}

inline double log_add(double a, double b)
{
	if(b>a)
	{
		double temp=a;
		a=b;
		b=temp;
	}
	return a+log(1+exp(b-a));
}

double norm(mat X)
{
	int n=X.n_rows;
	double sum=0;
	for(int i=0;i<n;i++)
	{
		sum+=abs<double>(X(i,0));
	}
	return sum;
}
double norm2(mat X)
{
	int n=X.n_rows;
	double sum=0;
	for(int i=0;i<n;i++)
	{
		sum+=pow(X(i,0),2);
	}
	return sum;
}

double mean2(mat X,double *w)
{
	int M=X.n_rows;
	double sum1=0;
	double sum2=0;
	for(int i=0;i<M;i++)
	{
		sum1+=X(i,0)*exp(w[i]);
		sum2+=exp(w[i]);
	}
	return sum1/sum2;
}

int is_present(mat X, mat Y)
{
	int res=0;
	int n=X.n_rows;
	for( int i=0;i<n;i++)
	{
		umat U=(X(i,span::all)==Y);
//		cout << U;
		if(as_scalar(sum(U,1))==U.n_elem)
		{
			res=1;
		}
	}
	return res;	
}

int nb_unique(mat X)
{
	int sum=0;
	int n=X.n_rows;
	mat Y=X(0,span::all);
	for(int i=1;i<n;i++)
	{
		if(!is_present(Y,X(i,span::all)))
		{
			sum++;
			Y=rbind<mat>(Y,X(i,span::all));			
		}
	}
return sum;	
}
mat Adjust(mat v, mat index)
{
	int p=v.n_rows;
	int nb=sum(sum(index));
	mat A(nb,nb);
	int k=0;
	A.fill(0);
	for(int i=0;i<p;i++)
	{
		if(index(i,0)==1)
		{
			A(k,k)=v(i,i);
			k++;
		}
	}
	return A;
}


int minvector(mat a)
{
	int n=a.n_elem;
	int k=0;
	double curr=a(0,0);
	for(int i=1;i<n;i++)
	{
		if(a(i,0)<curr)
		{
			k=i;
			curr=a(i,0);
		}
	}
	return k;
}

mat uglyGetridz(mat X)
{
	int p=X.n_cols;
	int n=X.n_rows;
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<p;j++)
		{
			if(abs(X(i,j))<0.0000001)	X(i,j)=0;
		}
	}
	return X;
}

mat del_mat(mat X, int k)
{
	int n=X.n_rows;
	int p=X.n_cols;
	mat Y(n-1,p-1);
	int m=0;
	int l=0;
	for(int i=0;i<n;i++)
	{
		if(i!=k)
		{
			l=0;
			for(int j=0;j<p;j++)
			{
				if(j!=k)
				{
					Y(m,l)=X(i,j);
					l++;
				}
			}
			m++;
		}
	}
	return Y;
}

mat del_vec(mat X, int k)
{
	int n=X.n_rows;
	mat Y(n-1,1);
	int m=0;
	for(int i=0;i<n;i++)
	{
		if(i!=k)
		{
			Y(m,0)=X(i,0);
			m++;
		}
	}
	return Y;
}
int Searchindex(imat x,int k)
{
	int l=-1;
	int _p=x.n_elem;
	for(int i=0;i<_p;i++)
	{
		if(x(0,i)==k) l=i;
	}
	return l;
}
#endif
